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solution.  In  this  research,  the  basic  theory  of  least-squares  finite  element  formulations  of  the  equations 
governing  viscous  incompressible  flows  and  shear  deformable  theories  of  plate  and  shell  structures  was 
carried  out  and  their  application  through  a  variety  of  benchmark  problems  was  illustrated.  In  the  case  of 
fluid  flows,  penalty  least-squares  finite  element  models  using  high  p-levels  and  low  penalty  parameters 
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TECHNICAL  REPORT 


Abstract 

This  report  summarizes  the  research  carried  out  under  Grant  F49620-03- 1-0201  on  the 
development  of  least-squares  based  finite  element  models  of  viscous  compressible  and 
incompressible  flows  as  well  as  shear  deformable  plates  and  shells.  The  main  objec¬ 
tive  of  this  research  was  to  develop  a  robust  and  accurate  computational  methodology 
based  on  least-squares  variational  principles  for  the  numerical  solution  of  the  equa¬ 
tions  governing  plates  and  shells  and  viscous  incompressible  and  compressible  fluid 
flows.  The  use  of  least-squares  principles  leads  to  a  variationally  unconstrained  min¬ 
imization  problem,  where  compatibility  conditions  between  approximation  spaces  - 
such  as  inf-sup  conditions  -  never  arise.  Furthermore,  the  resulting  linear  algebraic 
problem  will  always  have  a  symmetric  positive  definite  (SPD)  coefficient  matrix,  al¬ 
lowing  the  use  of  robust  and  fast  preconditioned  conjugate  gradient  methods  for  its 
solution.  In  this  research,  the  basic  theory  of  least-squares  finite  element  formulations 
of  the  equations  governing  viscous  incompressible  flows  and  shear  deformable  theories 
of  plate  and  shell  structures  was  carried  out  and  their  application  through  a  variety  of 
benchmark  problems  was  illustrated.  In  the  case  of  fluid  flows,  penalty  least-squares 
finite  element  models  using  high  p-lcvels  and  low  penalty  parameters  were  developed 
as  a  good  alternative  to  mixed  least-squares  finite  element  models,  also  developed  in 
the  research. 

The  new  computational  methodology  offers  many  theoretical  and  practical  advan¬ 
tages  over  traditional  finite  element  formulations.  In  the  context  of  solid  mechanics, 
least-squares  formulations  are  found  to  be  robust  for  the  bending  of  thin  and  thick 
plates,  effective  for  the  analysis  of  shell  structures  in  bending  and  membrane  domi¬ 
nated  states,  and  yield  accurate  predictions  of  generalized  displacements  and  stress 
resultants.  In  the  context  of  fluid  flows,  least-squares  formulations  (both  mixed  and 
penalty  models)  have  been  implemented  for  viscous  incompressible  fluid  flows.  Mixed 
least-squares  models  are  also  implemented  for  compressible  flows.  In  the  mixed  and 
penalty  finite  element  model  of  viscous  incompressible  flows,  high-order  expansions 
are  used  to  construct  the  discrete  model.  The  element-by-element  method  and  a 
matrix-free  version  of  the  conjugate  gradient  method  with  a  Jacobi  preconditioner 
are  used  to  solve  the  linear  system  of  equations.  Numerical  simulations  are  carried 
out  for  a  number  of  two-dimensional  benchmark  problems,  e.g.,  flow  over  a  backward 
facing  step  and  steady  flow  past  a  circular  cylinder.  The  effect  of  penalty  parameter 
on  accuracy  and  computational  cost  is  investigated  thoroughly  for  these  problems. 
The  least-squares  models  are  well  suited  for  large  scale  computations  and  robust  for 
moderately  high  Reynolds  number  and  high  Mach  number  flow  conditions  as  well  as 
for  thick  and  thin  structures. 
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We  find  that  the  A;- version  of  least-squares  (with  k  =  1),  achieves  equally  accurate 
results  when  compared  with  formulations  with  k  =  0,  at  a  lower  degree-of-freedom 
count.  However,  the  construction  of  the  k  =  1  basis  in  the  general  multi-dimensional 
setting  by  one-dimensional  tensor  products  has  undesirable  properties.  For  example, 
for  geometrically  distorted  elements  the  spectral  convergence  property  is  lost.  We 
proposed  in  this  research  the  weak-enforcement  of  k  =  1  continuity  through  the 
least-squares  functional.  In  other  words,  the  jump  of  the  primary  variables  and 
their  derivatives  across  inter-element  boundaries  is  minimized  in  a  least-squares  sense, 
through  the  least-squares  functional.  This  allows  the  use  of  practical  k  —  0  basis 
at  the  element  level,  and  we  achieve  k  —  1  continuity  globally.  In  addition,  the 
formulation  naturally  allows  for  geometric  and  basis  non-conformities  across  inter¬ 
element  boundaries. 


1  Introduction 

In  the  past  few  years  finite  element  models  based  on  least-squares  variational  princi¬ 
ples  have  drawn  considerable  attention  (see,  e.g.,  [1,2]).  In  particular,  given  a  partial 
differential  equation  (PDE)  or  a  set  of  partial  differential  equations,  the  least-squares 
method  allows  us  to  define  an  unconstrained  minimization  principle  so  that  a  finite 
element  model  can  be  developed  in  a  variational  setting.  The  idea  is  to  define  the 
least-squares  functional  as  the  sum  of  the  squares  of  the  equations  residuals  measured 
in  suitable  norms  of  Hilbert  spaces.  Assuming  the  governing  equations  (augmented 
with  suitable  boundary  conditions)  have  a  unique  solution,  the  least-squares  func¬ 
tional  will  have  a  unique  minimizer.  Thus,  by  construction,  the  least-squares  func¬ 
tional  is  always  positive  and  convex,  ensuring  coerciveness,  symmetry,  and  positive 
definitiveness  of  the  bilinear  form  in  the  corresponding  variational  problem.  More¬ 
over,  if  the  induced  energy  norm  is  equivalent  to  a  norm  of  a  suitable  Hilbert  space, 
optimal  properties  of  the  resulting  least-squares  formulation  can  be  established. 

However,  an  optimal  least-squares  formulation  may  result  in  an  impractical  finite 
element  model.  The  reconciliation  that  must  exist  between  practicality  and  opti¬ 
mality  in  least-squares  based  finite  element  models  is  of  great  importance  and  was 
first  recognized  by  Bochev  and  Gunzburger  [3,  2].  The  practicality  of  the  resulting 
finite  element  model  is,  to  a  large  extent,  determined  by  the  complexity  of  algorithm 
development  and  CPU  solve  time  of  the  resulting  discrete  system  of  equations.  Typ¬ 
ically,  the  practicality  is  measured  in  terms  of  Ck  continuity/regularity  of  the  finite 
element  spaces  across  inter-element  boundaries.  Ideally,  a  least-squares  finite  element 
model  with  “C°  practicality”  and  full  (mathematical)  optimality  is  to  be  developed 
-  unfortunately,  this  can  seldom  be  achieved  in  a  satisfactory  manner. 

Conforming  discretizations  require  that  the  finite  element  space  be  spanned  by 
functions  that  belong  to  the  Hilbert  space  H2m,  in  contrast  to  weak  form  Galerkin 
models  which  require  only  Hm  regularity  (due  to  the  weakened  differentiability  re¬ 
quirements  induced  by  the  integration  by  parts).  For  the  least-squares  model,  this 
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implies  a  minimum  of  C1  regularity  of  the  finite  element  spaces  across  inter-element 
boundaries  for  m  =  1. 

To  reduce  the  higher  regularity  requirements,  the  PDE  or  PDEs  are  first  trans¬ 
formed  into  an  equivalent  lower  order  system  by  introducing  additional  indepen¬ 
dent  variables,  sometimes  termed  auxiliary  variables,  and  then  formulating  the  least- 
squares  model  based  on  the  equivalent  lower  order  system.  The  additional  variables 
imply  an  increase  in  cost,  but  can  be  argued  to  be  beneficial  as  they  may  represent 
physically  meaningful  variables,  e.g.,  fluxes  or  stresses,  and  will  be  directly  approxi¬ 
mated  in  the  model.  Such  an  approach,  is  believed  to  be  first  explored  by  Jespersen  [4] 
and  is  the  preferred  approach  in  modern  implementations  of  least-squares  finite  el¬ 
ement  models.  For  2nd  order  PDEs,  an  equivalent  first-order  system  is  introduced, 
and  if  the  least-squares  functional  is  defined  in  terms  of  L 2  norms  only,  the  finite 
element  model  allows  the  use  of  nodal/modal  expansions  with  merely  C°  regularity. 

We  also  investigate  a  least-squares  formulation  where  the  “C°  practicality”  is 
relaxed  and  finite  element  spaces  are  allowed  to  retain  higher  regularity  across  inter¬ 
element  boundaries.  Such  formulations  do  not  require  that  auxiliary  variables  be 
introduced.  We  present  a  formulation  where  C 1  continuity  is  enforced  in  a  weak 
sense  through  the  least-squares  functional  thus  still  allowing  the  use  of  practical 
C°  nodal/modal  element  expansions.  The  latter  approach  naturally  allows  for  h- 
and  p-type  non-conformities  in  the  computational  domain  and  may  be  viewed  as  a 
discontinuous  least-squares  finite  element  formulation. 

In  this  report  we  give  a  brief  overview  of  the  work.  Specifically,  in  Section  2 
we  present  least-squares  formulations  using  a  unified  approach  for  an  abstract  initial 
boundary  value  problem,  which  could  represent  a  mathematical  model  describing 
fluid  flow  or  the  deformation  of  a  solid  structure.  In  Section  3  we  present  applications 
to  fluid  mechanics,  for  viscous  incompressible  flows,  and  in  Section  4  applications 
to  solid  mechanics,  for  shear-deformable  shell  structures.  We  refer  the  reader  to 
the  journal  publications  under  the  grant  for  additional  details  on  applications  to 
incompressible  and  compressible  flows  [5,  6,  7,  11,  12],  plates  and  shells  [8,  9],  and 
radiative  transfer  [10]. 

2  An  abstract  least-squares  formulation 

In  this  section  we  present  the  steps  involved  in  developing  and  arriving  at  a  least- 
squares  based  finite  element  model.  We  wish  to  present  the  procedure  in  a  general 
setting,  and  to  this  end  present  the  procedure  in  the  context  of  an  abstract  initial 
boundary  value  problem. 

2.1  Notation 

Let  O  be  the  closure  of  an  open  bounded  region  D  in  Rrf,  where  d  =  2  or  3  represents 
the  number  of  space  dimensions,  and  x  =  (aq, . . . ,  xd)  =  ( x ,  y,  z)  be  a  point  in  fi  = 
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U  dCl,  where  dCl  =  T  is  the  boundary  of  Q. 

For  s  >  0,  we  use  the  standard  notation  and  definition  for  the  Sobolev  spaces 
H*  (Q)  and  Hs  (r)  with  corresponding  inner  products  denoted  by  (•,  ■)SiSj  and  (•,  -)Sir 
and  norms  by  ||  •  ||Sin  and  ||  •  ||S) r,  respectively.  Whenever  there  is  no  chance  of 
ambiguity,  the  measures  Q  and  F  will  be  omitted  from  inner  product  and  norm 
designations.  We  denote  the  L2  (Q)  and  L2  (r)  inner  products  by  (•,  -)q  and  (•,  -)r, 
respectively.  By  Hs  (Q)  we  denote  the  product  space  [Hs  (fl)]d. 

2.2  The  abstract  problem 

Consider  the  following  abstract  initial  boundary  value  problem: 


A(u)  +  £x(u)  =  f 

infix  (0,  r] 

(1) 

G(u)  =  h 

on  T  x  (0,  t\ 

(2) 

in  which  Ct  and  £x  are  partial  differential  operators  in  time  and  space  respectively, 
acting  on  the  vector  u  of  unknowns.  For  example,  a  transient  scalar  Poisson  equation 
would  have  Ct(u)  =  du/dt  and  £x(«)  =  — V2u. 

The  vector  valued  function  f  is  a  known  forcing  function,  Q  is  a  trace  operator 
acting  on  u,  and  h  represents  a  known  vector  valued  function  on  the  boundary.  We 
assume  initial  conditions  are  given  such  that  the  problem  is  well  posed  and  a  unique 
solution  exists. 

2.3  A  first-order  system  least-squares  (FOSLS)  formulation 

The  L2  least-squares  functional  associated  with  the  abstract  initial  boundary  value 
problem  is  constructed  by  summing  up  the  squares  of  the  equations  residuals  in  the 
L-2  norm  and  is  given  by 

J  (u;f,h)  =  ^  (  ||  A(u)  +  £x(u)  -  f  ||oif2x(0,T]  +  II  0(u)  -  h  ||^rx(0iT] )  .  (3) 

It  is  easy  to  see  that  the  minimizer  of  (3)  solves  (l)-(2)  and  viceversa. 

Note  that  in  defining  the  FOSLS  functional  we  must  make  two  restrictions:  (1) 
the  temporal  and  spatial  partial  differential  operators  and  trace  operator  are  at  most 
of  first-order  and  (2)  the  least-squares  functional  is  defined  exclusively  in  terms  of  L2 
norms.  These  restrictions  are  necessary  in  order  to  ensure  a  pre-determined  level  of 
practicality  in  the  resulting  least-squares  based  finite  element  model:  specifically,  the 
permission  to  use  finite  element  spaces  with  merely  C°  regularity  across  inter-element 
boundaries. 

If  the  partial  differential  equations  (PDEs)  under  consideration  are  not  of  first- 
order,  the  “C'°  practicality”  of  the  least-squares  based  finite  element  model  comes  at 
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an  extra  cost,  implied  in  restriction  (1);  which  requires  that  the  partial  differential 
operators  be  of  first-order.  This  can  always  be  achieved  by  introducing  auxiliary 
variables  until  a  first-order  system  is  attained.  The  added  cost  might  be  viewed  as 
beneficial,  in  the  sense  that  the  auxiliary  variables  may  have  physical  relevance  to  the 
problem  under  consideration,  e.g.,  fluxes,  vorticity,  or  stresses. 

2.4  A  discontinuous  least-squares  (DLS)  formulation 

In  the  discontinuous  least-squares  formulation  the  functional  is  defined  at  the  element 
level,  so  that  jumps  across  inter-element  boundaries  may  be  minimized  in  the  least- 
squares  sense  and  global  weak  Ck,  k  >  0  is  achieved: 

We, 

J(u;f,h)  =  ^^e(ue;f,h); 

e—1 

Je  (ue; f.h)  =  |  (||  A(ue)  +  £x(ue)  -  f  Ho, oe x (o,t] 

+  ||  G( ue)  -  h  ||o,r=nrx(o,r]  +  II  ue  -  u°  Hl.r-x^.r]  )  (4) 

a 

where  the  superscript  a  denotes  abutting  elements  to  element  e. 

Unlike  the  FOSLS  formulation,  the  partial  differential  operators  need  not  be  re¬ 
duced  to  first-order  provided  the  index  k  in  the  last  residual  measure  minimizing 
jumps  across  inter-element  boundaries  is  at  minimum  k  =  n  —  1,  where  n  is  the  order 
of  the  partial  differential  operator.  In  addition,  practical  C°  basis  may  be  used  in  this 
formulation,  although  even  simpler  L2  subspaces  would  suffice.  For  n  >  1,  higher  con¬ 
dition  numbers  are  associated  with  this  formulation  due  to  the  higher  differentiability 
requirements. 

2.5  Time  stepping 

2.5.1  Space-time  coupled  approach 

Note  that  prior  to  defining  functionals  (3)  or  (4)  we  did  not  replace  the  temporal 
operator  with  a  discrete  equivalent.  This  results  in  a  fully  space-time  coupled  for¬ 
mulation,  implied  in  the  definition  of  functionals  (3)  or  (4)  where  the  L2  norm  is 
defined  in  space-time,  i.e.  ||  •  ||o,Ox(o,r]  denotes  the  L2  norm  of  the  enclosed  quantity 
in  space-time: 

T 

IMIo,nx(o,r]  ~  J  J  M  dt . 

o  n 

This  implies,  for  example,  that  a  two-dimensional  time-dependent  problem  will  be 
treated  as  a  three-dimensional  problem  in  space-time  domain.  When  dealing  with 
the  stationary  form  of  the  equations  the  integral  over  time  domain  is  simply  dropped. 
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In  the  space-time  coupled  approach,  the  effects  of  space  and  time  are  allowed  to 
remained  coupled.  There  is  no  approximation  of  the  initial  boundary  value  problem. 
Instead,  a  basis  is  introduced  in  time  domain  to  represent  the  time  evolution  of  the 
independent  variables. 

Invariably,  we  as  analysts  would  like  to  simulate  and  study  the  time  evolution  of 
an  initial  boundary  value  problem  for  large  values  of  time.  Taking  into  considera¬ 
tion  modelling  issues,  we  realize  that  this  would  require  a  space-time  mesh  with  a 
large  number  of  elements  in  time  domain.  The  size  of  the  resulting  set  of  assembled 
algebraic  equations  could  be  large  and  prohibitively  expensive  in  terms  of  available 
computer  memory  and  non-optimal  in  terms  of  CPU  solve  time.  To  alleviate  the 
drawbacks,  we  adopt  a  time-stepping  procedure  in  which  the  solution  is  obtained  for 
space-time  strips  in  a  sequential  manner.  The  initial  conditions  for  the  current  space- 
time  strip  are  obtained  from  the  latest  space  plane  from  the  previous  space-time  strip. 
Hence,  for  each  space-time  strip  we  solve  a  true  initial  boundary  value  problem,  by 
minimizing  the  following  (e.g.  FOSLS)  functional  in  space-time  domain: 

J(uif,h)=i(||£,(u)  +  £x(u)-f||'„>,Miti| 

+  ll£(u)-h|lo1rxM»+i])  (5) 

where  the  interval  [tS)  t3+i]  can  be  taken  arbitrarily  large,  i.e.,  there  are  no  restrictions 
on  the  size  of  the  interval.  Additional  details  on  this  time  stepping  approach  may  be 
found  in  Pontaza  and  Reddy  [6]. 

2.5.2  Space-time  decoupled  approach 

Alternatively,  the  temporal  operator  can  be  represented  by  truncated  Taylor  series 
expansions  in  time  domain,  e.g.  a  backward  Euler  or  trapezoidal  rule  approximation. 
First,  the  temporal  operator  in  Eq.  (1)  is  replaced  by  the  discrete  approximation: 

A ( u)  «  £A<(us+1,  us_<2) ,  q  =  0, 1, 2, . . . 

where  the  time  increment  dependence  of  the  discrete  operator  is  explicit  as  well  as 
its  dependence  on  histories  of  previous  time  steps.  For  example,  for  Ct(u)  —  du/dt  a 
(backward  Euler)  first-order  approximation  would  be  £At(u)  —  (us+1  —  us)/At. 

To  march  the  problem  in  time,  we  must  minimize  the  following  (e.g.  FOSLS) 
functional  at  each  time  step: 

Ju  (u;  f,  h)  =  i  (  ||  £At( iT+1,  iF-*)  +  £x(us+1)  -  fs+l  ||g.n 

+  II  £(us+1)  —  hs+1  |lo,r  )  (6) 

where  the  dependence  on  the  time  increment  At  =  ts+ j  —  ts  is  evident.  This  time 
stepping  approach  obviously  has  associated  with  it  a  much  lower  computational  cost 
when  compared  with  the  space-time  coupled  approach,  as  the  dimensionality  of  the 
problem  is  not  increased. 
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2.6  The  variational  problem 

Having  defined  the  least-squares  functional,  the  abstract  least-squares  minimization 
principle  can  be  stated  as: 

find  u£X  such  that  J  (u;  f,  h)  <  J  (v;  f,h)Vv£X  (7) 

where  X  is  a  suitable  vector  space,  e.g.  X  —  H1  (Q)  for  a  FOSLS  space-time  decou¬ 
pled  formulation. 

The  Euler-Lagrange  equation  for  this  minimization  problem  is  given  by  the  fol¬ 
lowing  variational  problem: 

find  u£l  such  that  B  (u,  v)  =  T  (v)  V  v  £  X  (8) 

where  B  is  a  symmetric  form  given  by 

H(u,v)  =  (£(u),£(v))n  +  (£(u),  £(v))r 
and  J7  is  a  functional  given  by 

^(v)  =  (f,£(v))n  +  (h,£(v))r 


where  C  —  £At  +  £x- 

The  inclusion  of  the  boundary  residual  in  the  least-squares  functional  allows  the 
use  of  spaces  X  that  are  not  constrained  to  satisfy  the  boundary  condition  (2).  In 
such  a  case,  the  boundary  condition  (2)  is  enforced  in  a  weak  sense  through  the  least- 
squares  functional.  This  is  a  tremendous  advantage  of  least-squares  based  formula¬ 
tions,  as  it  allows  boundary  conditions  that  are  computationally  difficult  to  impose  to 
be  efficiently  included  in  the  least-squares  functional.  An  example  where  this  prop¬ 
erty  becomes  extremely  useful  is  for  viscous  or  inviscid  compressible  flow,  where 
characteristic-based  boundary  conditions  need  to  be  prescribed  at  outflow/inflow 
boundaries.  Of  course,  if  the  boundary  condition  (2)  can  be  easily  imposed  and 
included  in  the  space  X,  we  omit  the  residual  associated  with  the  boundary  term  in 
the  least-squares  functional. 

The  abstract  expressions  given  above  for  the  symmetric  form  B  and  functional  T 
are  only  valid  when  the  partial  differential  operators  are  linear.  For  the  case  when 
the  partial  differential  operators  are  nonlinear,  the  following  more  general  expressions 
apply 

B{  u,v)  =  (£(u),6£(u))n  +  (£(u),  6£(u))r 

and 

JF(v)  =  (f,5£(u))n  +  (h,5C?(u))r 

where  it  is  understood  that  <5  u  =  v.  In  general,  when  the  partial  differential  operators 
are  nonlinear,  the  resulting  form  will  be  non-symmetric.  Symmetry  of  the  form  is 
restored  only  when  the  Euler-Lagrange  equation  is  linearized  by  Newton’s  method. 
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2.7  The  finite  element  model 


The  finite  element  model  is  obtained  by  either  restricting  (8)  to  the  finite  dimensional 
subspace  Xhp  of  the  infinite  dimensional  space  X,  or  equivalently  by  minimizing  (3) 
with  respect  to  the  chosen  approximating  spaces.  This  process  leads  to  the  discrete 
variational  problem  given  by 

find  uhp  €  Xhp  such  that  B  (uhp,  vhp)  =  T  (vfcp)  V  vhp  G  Xhp  (9) 

We  proceed  to  define  a  discrete  problem  by  choosing  appropriate  finite  element 
subspaces  for  each  of  the  components  of  the  vector  valued  function  u.  There  are 
no  restrictive  compatibility  conditions  on  the  discrete  spaces,  so  we  choose  the  same 
finite  element  subspace  for  each  of  the  primary  variables. 


2.8  Nodal/modal  expansions 


We  present  in  this  section  some  details  on  the  high-order  nodal/modal  expansions 
used  in  this  work.  Advantages  of  using  high-order  methods  include  [13,  5,  6]:  ex¬ 
ponentially  fast  decay  of  error  measures  for  smooth  solutions,  small  diffusion  and 
dispersion  errors,  better  data  volume  over  surface  ratio  allowing  for  high  efficien¬ 
cies  in  parallel  processing,  and  higher  efficiency/accuracy  in  long  time  integration  of 
unsteady  problems. 

In  a  modal  expansion,  the  finite  element  spaces  are  spanned  by  tensor  products 
of  the  one-dimensional  C°  p-type  hierarchical  basis 


i’iiO  = 


'Izi 

<  (¥)  (¥)  ptf 

i+i 

2 


1  =  1 

2  <i  <p,  p  >  2 

i  =  p+  1 


(10) 


In  definition  (10),  Pp'^  are  the  Jacobi  polynomials  of  order  p.  We  use  ultraspheric 
polynomials  corresponding  to  the  choice  a  —  0  with  a  =  (3  =  0  or  1. 

Figure  1  shows  the  one-dimensional  modal  basis  for  the  case  p  =  5.  The  linear 
basis  or  “hat-functions”  ensure  the  C°  continuity  requirement  across  element  bound¬ 
aries  and  the  p-bubbles  hierarchically  enrich  the  finite  element  space.  Note  that  by 
construction  the  p-bubbles  vanish  at£  =  — 1,£  =  +1  and  have  no  nodes  associated 
with  them. 


In  a  nodal  expansion,  the  finite  element  spaces  are  spanned  by  tensor  products  of 
the  one-dimensional  C°  spectral  nodal  basis 


(^-i)(^  +  i)l;(0 

p(p  +  l)Lp(&)(£-&) 


(11) 


In  Eq.  (11),  Lp  =  Pp’°  is  the  Legendre  polynomial  of  order  p  and  &  denotes  the 
location  of  the  roots  of  (£  —  1)(£  +  !)!£(£)  =  0  in  the  interval  [—1, 1].  The  set  of 
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Figure  1:  C°  p-type  hierarchical  modal  basis.  Shown  is  the  case  of  p  —  5.  The 
p-bubbles  are  scaled  by  a  factor  of  4,  for  viewing  ease. 


Figure  2:  C°  p- type  (spectral)  nodal  basis.  Shown  is  the  case  of  p  =  4. 
points  are  commonly  referred  to  as  the  Gauss-Lobatto-Legendre  (GLL)  points. 

Figure  2  shows  the  one-dimensional  nodal  basis  for  the  case  p  =  4.  The  location 
of  the  nodes  coincides  with  the  roots  of  the  aforementioned  Legendre  polynomial  and 
thus  receives  the  name  of  a  “spectral”  basis.  The  Kronecker  delta  property  is  evident 
from  the  figure  and  is  an  attractive  feature  of  this  basis,  as  the  coefficients  coincide 
with  nodal  values. 


3  Viscous  incompressible  fluid  flows 

The  numerical  solution  of  the  incompressible  Navier-Stokes  equations  using  least- 
squares  based  finite  element  models  is  among  the  most  popular  applications  of  least- 
squares  methods.  Least-squares  formulations  for  incompressible  flow  circumvent  the 
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inf-sup  condition,  thus  allowing  equal-order  interpolation  of  velocities  and  pressure, 
and  result  (after  suitable  linearization)  in  linear  algebraic  systems  with  a  SPD  coef¬ 
ficient  matrix.  This  translates  into  easy  algorithm  development  and  leads  to  the  use 
of  robust  and  fast  iterative  solvers,  resulting  in  substantial  improvements  over  the 
traditional  weak  form  Galerkin  finite  element  model  -  where  the  finite  element  spaces 
for  velocities  and  pressure  must  satisfy  an  inf-sup  compatibility  condition  and  one 
must  deal  with  an  un-symmetric  and  indefinite  coefficient  matrix. 

3.1  The  incompressible  Navier- Stokes  equations 

We  consider  the  solution  of  the  Navier-Stokes  equations  governing  incompressible 
flow,  which  in  dimensionless  form  can  be  stated  as  follows: 

Find  the  velocity  u  (x,  t)  and  pressure  p  (x,  t)  such  that 


(?U  .  1 

^  +  (u.V)u+VP--Vu  =  f 

in  fl  x  (0,  r] 

(12) 

<1 

C 

II 

o 

in  fl  x  (0,  t) 

(13) 

u  (x,  0)  =  °u  (x) 

in  fl 

(14) 

u  =  us 

on  r„  x  (0,  r] 

(15) 

n  •  o-  -  f 5 

onT/  x  (0,  r] 

(16) 

where  T  =  U  T/  and  Tu  fl  T/  =  0.  Re  is  the  Reynolds  number,  V  ■  °u  =  0,  f 
is  a  dimensionless  force,  n  is  the  outward  unit  normal  on  the  boundary  of  fl,  lT 
is  the  prescribed  velocity  on  the  boundary  r„,  fs  are  the  prescribed  tractions  on 
the  boundary  Tf ,  and  in  Eq.  (14)  the  initial  conditions  are  given.  The  conditions 
on  the  boundary  Tf  in  Eq.  (16)  are  used  to  model  outflow  conditions,  with  a  = 
—pi  +  (1/Re)  Vu  and  fs  =  0. 

As  discussed  in  the  previous  section,  in  a  FOSLS  formulation  the  governing  equa¬ 
tions  must  be  recast  as  an  equivalent  first-order  system.  This  will  allow  the  use  of 
practical  C°  basis  in  the  finite  element  model.  We  use  a  vorticity-based  first-order 
system,  using  u  —  V  x  u,  and  in  view  of  the  vector  identity 

V  x  V  x  u  =  —  V2u  +  V  (V  •  u) 

and  the  incompressibility  constraint  given  in  Eq.  (13),  the  non-stationary  Navier- 
Stokes  equations,  Eqs.  (12)-(16),  can  be  replaced  by  their  first-order  system  equiva¬ 
lent: 

Find  the  velocity  u  (x,  t),  pressure  p  (x,  t),  and  vorticity  u)  (x,  t )  such  that 

X 

-^7  +  (u  •  V)  u  +  Vp  +  —  V  x  u  =  f  inflx(0,T]  (17) 

ot  Ke 

w-Vxu  =  0  in  fl  x  (0,  r]  (18) 
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V  •  u  =  0 

in  D  x  (0,  r] 

(19) 

o 

il 

3 

t> 

in  Q  x  (0,  r] 

(20) 

u  (x,  0)  =  °u  (x) 

in  Q 

(21) 

u  =  us 

on  Tu  x  (0,  t] 

(22) 

u  =  us 

on  x  (0,  r] 

(23) 

h  ■  cr  =  fs 

on  Tf  x  (0,  r] 

(24) 

where  D  rw  =  0,  i.e.  if  velocity  is  specified  at  a  boundary,  vorticity  need  not  be 
specified  there.  This  implies  that  no  artificial  boundary  conditions  for  vorticity  need 
to  be  devised  at  boundaries  where  the  velocity  is  specified. 

Other  first-order  systems  are  also  possible,  e.g.  a  stress-based  first-order  system 
or  a  velocity  gradient  based  first-order  system.  Recasting  the  governing  equations  as 
a  first-order  system  is  not  necessary  when  using  the  DLS  formulation. 


3.2  Numerical  examples 

3.2.1  Kovasznay  flow 


The  first  benchmark  problem  to  be  used  for  the  purposes  of  verification  is  an  ana¬ 
lytic  solution  to  the  two-dimensional,  stationary  incompressible  Navier-Stokes  due  to 
Kovasznay  [14].  The  spatial  domain  in  which  Kovasznay’s  solution  is  defined  is  taken 
here  as  the  bi-unit  square  Q  =  [—0.5, 1.5]  x  [—0.5, 1.5].  The  solution  is  given  by 

u(x,  y)  =  1  —  eXx  cos(27 ry) 


v(x,y)  =  —  eXx  sin(27r y)  (25) 

P(®,  V)  =  Po  ~  \  e2Ax 


where  A  =  Re/2  —  (Re2/4  +  47r2)1/2,  p0  is  a  reference  pressure  (an  arbitrary  constant), 
and  we  choose  Re  =  40. 

Figure  3  shows  p-convergence  curves  of  the  velocity  field  in  the  H1  norm  as  a 
function  of  total  number  of  degrees  of  freedom  for  the  different  formulations. 

We  see  that  spectral  convergence  of  the  velocity  field  is  realized  for  all  the  formu¬ 
lations.  However,  the  FOSLS  formulations  have  a  higher  degree  of  freedom  count  due 
to  the  auxiliary  variables  used  to  recast  the  governing  equations  as  an  equivalent  first- 
order  system.  From  the  numerical  results  it  might  appear  that  the  DLS  formulation  is 
the  formulation  of  choice.  However,  due  to  the  higher  order  operators  involved  in  the 
DLS  formulation  (e.g.  the  V2  operator),  the  conditioning  of  the  resulting  coefficient 
matrix  is  higher. 
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Figure  3:  Convergence  of  the  velocity  field  to  the  Kovasznay  solution  in  the  H 1  norm 
for  FOSLS  and  DLS  formulations. 

3.2.2  Flow  past  two  circular  cylinders  in  a  side-by-side  arrangement 

We  consider  two-dimensional  flow  of  an  incompressible  fluid  past  two  circular  cylinders 
in  a  side-by-side  arrangement.  Both  cylinders  are  equal  in  size,  with  diameter  D,  and 
face  the  free-stream.  The  flow  around  such  an  arrangement  is  characterized  by  three 
distinct  flow  regimes,  depending  on  the  gap  size  S  between  cylinder  surfaces  [15,  16]. 

The  cylinders  are  of  unit  diameter  and  are  at  a  distance  S/D  —  0.85  of  each  other. 
The  simulation  is  carried  out  using  a  space-time  coupled  formulation.  The  connected 
model  in  space-time,  Clh  x  [ts,  ta+1],  consists  of  762  finite  elements  in  space  and  a  single 
element  layer  in  time.  Figure  4  shows  the  connected  model  in  space  and  a  close-up 
view  of  the  geometric  discretization  around  the  circular  cylinders. 

We  use  nodal  expansions  with  =  4  and  p7  =  2  in  each  element,  i.e. 

fourth-order  expansions  in  space  and  a  second-order  expansion  in  time,  resulting  in 
Wlof  =  149, 508  for  a  space-time  strip.  At  each  Newton  step  the  linear  system  of 
algebraic  equations  is  solved  using  the  matrix-free  conjugate  gradient  algorithm  with 
a  Jacobi  preconditioner.  For  the  time  marching  procedure  the  size  of  the  time  step, 
At  =  ts+ 1  —  ts,  was  chosen  as  At  =  0.20.  We  consider  a  Reynolds  number  of  100, 
based  on  the  free-stream  velocity  and  cylinder  diameter. 

At  the  upstream  boundary  of  the  computational  domain  both  velocity  components 
are  assigned  free-stream  values:  u  =  Uoo  =  1  and  v  —  0.  At  the  lateral  boundaries 
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(a)  (b) 


x 


Figure  4:  Computational  domain  and  mesh  for  flow  past  two  circular  cylinders  in  a 
side-by-side  arrangement,  S/D  =  0.85.  (a)  Connected  model,  Clh.  (b)  Close-up  view 
of  the  geometric  discretization  around  the  circular  cylinders. 


a  no-flux  boundary  condition  is  imposed:  du/dy  —  0  and  v  =  0.  No-slip  boundary 
conditions  are  specified  at  the  cylinder  surface:  u  —  v  =  0.  The  outflow  boundary 
condition  is  imposed  in  a  weak  sense  through  the  least-squares  functional. 

At  this  gap  size,  we  did  not  see  (or  expected)  a  well-defined  periodic  steady  state. 
The  simulation  was  thus  carried  out  for  t  €  [0,400],  by  which  time  the  flow  exhibited 
“well-developed”  characteristics  such  as  intermittent  bistable  gap  jets  and  amalga¬ 
mation  of  gap  vortices  leading  to  the  formation  of  a  single  large-scale  vortex  street. 

Figure  5  shows  instantaneous  vorticity  contours  at  (a)  t  =  314  and  (b)  t  —  352, 
at  which  times  the  gap  jet  is  “shooting”  upwards  and  downwards  respectively.  In 
accordance  with  the  experimental  visualizations  of  Williamson  [15],  a  single  large-scale 
vortex  street  is  formed  downstream  of  the  cylinders,  by  virtue  of  vortex  interactions 
in  the  near-wake  of  the  cylinders.  Gap  vortices  from  both  cylinders  are  squeezed  and 
amalgamated  with  dominant  outer  vortices,  predominantly  towards  the  narrow-wake 
side  (the  gap  jet  “shooting”  direction). 

In  Fig.  6  we  plot  the  force  coefficient  associated  with  the  repulsive  force  expe¬ 
rienced  by  the  circular  cylinder  whose  center  is  located  at  ( x,y )  =  (0,-0.925).  At 
early  times,  0  <  t  <  125,  we  see  a  well-defined  shedding  frequency,  due  to  shedding 
synchronization  in  antiphase.  For  t  >  150,  the  flow  becomes  asymmetric  due  to  the 
bistable  biased  gap  flow. 

The  value  of  the  L2  least-squares  functional  for  remained  below  10~3  throughout 
the  time  marching  procedure,  meaning  that  conservation  of  mass  and  momentum 
are  being  satisfied  to  within  10~3  at  all  times  -  implying  a  time-accurate  numerical 
simulation. 
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t  =  314.00 


Figure  5:  Instantaneous  vorticity  contours  for  the  flow  around  two  circular  cylinders 
in  a  side-by-side  arrangement  with  a  gap  size  of  S/D  =  0.85  and  Re^  =  100. 


Figure  6:  Time  history  of  repulsive  (lift)  coefficient  experienced  by  the  circular  cylin¬ 
der  whose  center  is  located  at  (x,y)  —  (0,  —0.925). 
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4  Shear-deformable  shells 


Finite  element  formulations  for  the  analysis  of  plates  and  shell  structures  are  tradition¬ 
ally  derived  from  the  principle  of  virtual  displacements  or  the  principle  of  minimum 
total  potential  energy  [17,  18].  When  considering  the  limiting  behavior  of  a  shell  as 
the  thickness  becomes  small,  for  a  given  shell  geometry  and  boundary  conditions,  the 
shell  problem  will  in  general  fall  into  either  a  membrane  dominated  or  bending  dom¬ 
inated  state  -  depending  on  whether  the  membrane  or  bending  energy  component 
dominates  the  total  energy.  Displacement-based  finite  element  models  have  no  major 
difficulties  in  predicting  the  asymptotic  behavior  of  the  shell  structure  in  the  mem¬ 
brane  dominated  case.  However,  computational  difficulties  arise  in  the  case  when  the 
deformation  is  bending  dominated.  A  strong  stiffening  of  the  element  matrices  occurs, 
resulting  in  spurious  predictions  for  the  membrane  energy  component.  This  phenom¬ 
enon  is  known  as  membrane-locking.  In  shear-deformable  shell  models,  yet  another 
form  of  locking  occurs  and  presents  itself  (again)  in  a  strong  stiffening  of  the  element 
matrices,  resulting  in  spurious  predictions  for  the  shear  energy  component.  This  form 
of  locking  is  also  present  in  plate  bending  analysis  when  the  side-to-thickness  ratio 
of  the  plate  is  large  (i.e.,  when  modelling  thin  plates).  This  locking  phenomenon  is 
known  as  shear-locking. 

Least-squares  finite  element  formulations  for  plates  and  shells  have  been  shown 
to  be  robust  with  regards  to  membrane-  and  shear-locking  and  to  yield  highly  ac¬ 
curate  results  for  displacements  as  well  as  stresses  (or  stress  resultants)  [8,  9].  The 
formulations  retain  the  generalized  displacements  and  stress  resultants  as  indepen¬ 
dent  variables  and,  in  view  of  the  nature  of  the  variational  setting  upon  which  the 
finite  element  model  is  built,  allows  for  equal-order  interpolation.  In  the  following 
we  present  numerical  results  for  a  barrel  vault  (cylindrical  shell)  loaded  by  its  own 
weight. 


4.1  Governing  equations 

We  consider  circular  cylindrical  shells,  where  the  shell  mid-surface  S  is  given  by 

S  —  {—L  <  xi  <  L ,  x\  +  £3  =  R2  |  (xi,  x2,x3)  £  R3}  C  R3  ,  (26) 

where  2 L  and  R  are  the  length  and  radius  of  the  shell.  The  shell  mid-surface  S,  given 
by  Eq.  (26),  can  be  parametrized  by  the  single  chart  (j)  =  </>2,  03 ),</>:  Q  C  R2  — > 

S  C  IP:3, 

<M£\£2) 

=R  M?/R)  (27) 

m\e)  =r  cos  (e/R) 

so  that  is  the  rectangle  occupying  the  region 

{(£\£2)  G  Q  |  -  L  <  ?  <  L,  Rtt  <  i2  <  Rtt}  C  R2  .  (28) 
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In  Naghdi’s  shear-deformable  shell  model  [19,  20],  the  membrane,  bending,  and 
shear  strain  measures  (eap,  Xa[h  («)  are  [9] 

£ll  =  uhl  ,  2  £12  =  1^1,2  +  U2<1  ,  £22  =  ^2,2  +  (29) 

Xll  =  01,1  ,  2  Xi2  =  01,2  +  $2,1  +  ,  X22  —  02,2  +  —  (^2,2  +  (30) 


Cl  —  U3,l  +  01  >  C2  —  ^3,2  +  02  — 


and  the  equilibrium  equations  take  the  form 
iV1/  +  JVf  +  p1  =  0 


„  99  M\2  Mf  Q2  . 

N?  +  Nf +  +  ^+p2  =  o 


R  +  R 


Q\  +  Q\ 2 


AT 
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M22 


+  p3  =  0 


(31) 

(32) 

(33) 

(34) 

(35) 


(36) 


R  R2 

M\l  +  Mf  -  Q1  =  0 
Mj  +  Mf  -  Q2  =  0 

where  ua  are  the  displacements  of  the  shell  mid-surface,  113  is  the  out-of-plane  dis¬ 
placement,  0 a  are  rotations  of  the  transverse  material  fibers  originally  normal  to  the 
shell  mid-surface,  and  Na3,  Ma3,  Q°  are  the  membrane,  bending,  and  shear  thickness- 
averaged  stress  resultants.  Here  we  employ  the  convention  that  Greek  indices  range 
over  1  and  2  and  that  denotes  differentiation. 

The  stress  resultants  are  related  to  the  strain  measures  through  the  following 
constitutive  relations  [19] 


N 


11 


tE  .  _r22  tE  . 

(£11  +  ^£22) ,  N  —  7; - 777  (^  £11  +  £22) 


(1-^ 


(1  -  S) 


N 12  = 


tE 


(1  +  *) 


£12 


(37) 


M11  = 


t3  E 

12  (1  —  v2) 


(xii  T  1^X22)  j  M 
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t3E 

12(1  -  v 2) 


{vXn  +  X22) 


M12  = 


t3E 


12(1  +  1/) 


X12 


Q1 


tE 


Ks  Ci,  Q2 


A's  (2 


(38) 


(39) 


2(l  +  i/)“s"1’  2(1  +  */) 

where  t  is  the  thickness  of  the  shell,  E  is  the  Young’s  modulus,  v  is  the  Poisson’s 
ratio,  and  Ks  is  the  shear  correction  factor  for  the  isotropic  material.  If  we  let  R  — »  00 
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we  recover  the  (linear)  shear-deformable  plate  bending  strain  measures  and  governing 
equations,  where  membrane  and  bending  effects  are  decoupled. 

The  equilibrium  equations  (32)-(36)  and  constitutive  relations  (37)-(39)  are  al¬ 
ready  of  first-order  and  are  used  to  define  the  least  squares  functional.  The  least- 
squares  formulation  and  finite  element  model  follow  from  the  outline  given  in  Sec¬ 
tion  2. 


4.2  Numerical  example:  Barrel  vault 

We  consider  a  barrel  vault  loaded  by  its  own  weight.  The  barrel  vault  is  a  segment 
of  a  circular  cylindrical  shell  whose  mid-surface,  after  being  parametrized  by  (27),  is 
given  by 

n={(e,e)\-L<e<L,-^R<e<2-iiR}.  (40) 

The  barrel  vault  is  simply-supported  on  rigid  diaphragms  on  opposite  edges  and  is 
free  on  the  other  two  edges.  For  the  described  loading,  geometry,  and  boundary 
conditions,  the  problem  is  popularly  known  as  the  Scordelis-Lo  roof. 

By  symmetry  considerations,  the  computational  domain  is  limited  to  1/4  of  the 
total  shell,  so  that 


nk  =  {(e,e)\o<e  <L,o<e  <2i-R}-  m 


The  geometry  of  the  barrel  vault  is  specified  as  follows:  2 L  =  50  ft,  R  =  25  ft, 
and  t  =  3  in.,  so  that  R/t  —  100.  The  material  is  homogeneous  and  isotropic  with 
E  —  3  x  106  psi  and  u  =  0.  The  shear  correction  factor  Ks  is  specified  as  5/6  and 
the  self-weight  loading  as  pz  =  90  lb/ft2  uniformly  distributed  over  the  surface  area 
of  the  vault. 

The  connected  model,  flk  C  R2,  consists  of  a  4  x  4  finite  element  mesh.  For 
illustrative  purposes  we  present  in  Fig.  7  the  finite  element  mesh  on  the  entire  mid¬ 
surface  of  the  barrel  vault,  <S  C  R3.  The  mesh  is  regular  (i.e.,  not  distorted)  and 
graded.  We  expect  strong  boundary  layers  in  the  stress  resultant  profiles  along  the 
free  and  supported  edges,  so  the  mesh  is  graded  towards  those  regions. 

First,  we  present  a  convergence  study  in  strain  energy  for  increasing  p- levels  of 
the  element  approximation  functions.  An  analytic  value  for  the  strain  energy  is  not 
available,  so  we  use  instead  a  reference  value.  The  reference  value  was  obtained  using 
displacement  based  weak  form  Galerkin  elements  with  a  p- level  of  12.  Denoting  by 
Ure{  the  reference  strain  energy  of  the  barrel  vault,  the  error  measure  is  given  by 


E  = 


|  ^ref  _  yhp 

U™{ 


(42) 


In  Fig.  8  we  plot  the  error  measure  E  for  the  least-squares  formulation  as  a  function 
of  the  expansion  order  in  a  logarithmic-linear  scale.  We  see  from  the  convergence  in 
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strain  energy  curve  that  an  accurate  least-squares  solution  is  achieved  for  p-levels  of  6 
and  higher. 

Table  1  shows  results  for  the  vertical  displacement  and  stress  resultants  at  the 
center  of  the  free  edge  of  the  barrel  vault,  (£\£2)  =  (0,  qf-ft),  for  p-levels  of  4,  6,  8, 
and  10.  Similarly,  in  Table  2  we  present  results  for  the  vertical  displacement  and  stress 
resultants  at  the  crown  of  the  barrel  vault,  (£*,  £2)  =  (0, 0).  We  see  from  the  tabulated 
data  that  the  change  in  point  values  at  p-levels  of  6,  8,  and  10  is  negligible,  and  thus 
a  converged  numerical  solution  could  be  declared  at  p- levels  of  6  or  8.  The  predicted 
vertical  deflection  at  the  center  of  the  free  edge  (Table  1)  is  in  good  agreement  with 
the  shallow  shell  analytical  value  of  3.7032  in.,  the  commonly  used  reference  value  for 
finite  element  analysis  of  3.6288  in.  [21],  and  the  value  of  3.6144  in.  obtained  using 
an  assumed  strain  method  in  a  fine  mesh  [22]. 


z 


Figure  7:  Finite  element  mesh  on  the  entire  mid-surface  of  the  barrel  vault,  S  C  K3, 
showing  the  surface  coordinate  system  (£1,£2)  6  M2. 
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Figure  8:  Convergence  of  strain  energy  for  the  barrel  vault  problem. 


Table  1:  p-convergence  study  showing  vertical  displacement  and  stress  resultants  at 
the  center  of  the  free  edge  of  the  barrel  vault. 


p  level 

w  (in.) 

Nn  (kip/ft) 

M11  (kip  ft/ft) 

4 

-3.1208 

68.3942 

-0.5610 

6 

-3.6162 

75.7476 

-0.6400 

8 

-3.6173 

75.7582 

-0.6400 

10 

-3.6174 

75.7593 

-0.6400 

Table  2:  p-convergence  study  showing  vertical  displacement  and  stress  resultants  at 
the  crown  of  the  barrel  vault. 


p  level 

w  (in.) 

TV11  (kip/ft) 

N22  (kip/ft) 

Mn  (kip  ft/ft) 

M22  (kip  ft/ft) 

4 

-3.4148 

0.0714 

1.7597 

6 

-1.5835 

-3.4861 

0.0959 

8 

-3.4862 

0.0959 

10 

-3.4862 

0.0959 
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5  Scientific  Progress  and  Accomplishments 

A  new  computational  methodology  based  on  least-squares  variational  principles  and 
the  finite  element  method  is  developed  for  the  numerical  solution  of  the  stationary  and 
non-stationary  Navier-Stokes  equations  governing  viscous  incompressible  and  com¬ 
pressible  fluid  flows  and  nonlinear  equations  governing  shear  deformation  theories  of 
plate  and  shell  structures.  The  use  of  least-squares  principles  leads  to  a  variationally 
unconstrained  minimization  problem,  where  compatibility  conditions  between  approx¬ 
imation  spaces  -  such  as  inf-sup  conditions  -  never  arise.  Furthermore,  the  resulting 
linear  algebraic  problem  will  always  have  a  symmetric  positive  definite  (SPD)  coef¬ 
ficient  matrix,  allowing  the  use  of  robust  and  fast  preconditioned  conjugate  gradient 
methods  for  its  solution.  In  the  context  of  viscous  incompressible  flows,  least-squares 
based  formulations  offer  substantial  improvements  over  the  (traditional)  weak  form 
Galerkin  finite  element  models  -  where  the  finite  element  spaces  for  velocities  and 
pressure  must  satisfy  an  inf-sup  compatibility  condition  and  one  must  deal  with  an 
unsymmetric  and  indefinite  coefficient  matrix.  In  contrast,  least-squares  formulations 
circumvent  the  inf-sup  condition,  thus  allowing  equal-order  interpolation  of  velocities 
and  pressure,  and  result  (after  suitable  linearization)  in  linear  algebraic  systems  with 
a  SPD  coefficient  matrix. 

A  penalty  least-squares  finite  element  model  is  also  developed  as  a  better  alter¬ 
native  to  traditional  penalty  finite  element  model.  Advantage  of  the  penalty  least- 
squares  finite  element  model  is  that  it  gives  very  accurate  results  for  very  low  penalty 
parameters  when  used  with  high  order  element  expansions.  It  is  found  that  the  com¬ 
puted  pressure  fields  are  continuous,  and  their  values  are  found  to  be  in  excellent 
agreement  with  published  results. 

We  have  also  developed  a  least-squares  formulation,  where  regularity  of  order  k 
is  achieved  by  the  weak  enforcement  of  continuity  constraints  across  inter-element 
boundaries.  This  allows  for  the  use  of  practical  k  —  0  expansions  at  the  element 
level,  and  can  achieve  any  desired  regularity  at  the  global  level.  The  formulation 
naturally  allows  for  h—  and  p— type  non-conformities. 

The  following  research  has  been  accomplished: 

•  Developed  mixed  least-squares  finite  element  models  of  the  Navier-Stokes  equa¬ 
tions  governing  viscous  incompressible  flows. 

•  Developed  space-time  coupled  least-squares  finite  element  models  of  non-stationary 
Navier-Stokes  equations  governing  viscous  incompressible  flows. 

•  Developed  least-squares  finite  element  models  of  equations  governing  viscous 
compressible  flows. 

•  Developed  least-squares  finite  element  models  of  bending  of  shear  deformable 
plates  and  shells. 
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•  Developed  penalty  least-squares  finite  element  models  of  the  stationary  Navier- 
Stokes  equations  governing  viscous  incompressible  flows. 

•  Developed  weak  k-ve rsion  least-squares  finite  element  models,  allowing  for  h- 
and  p-type  nonconformity. 
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